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Express Mail #: EL069309088US 



ADAPTIVE TRAIN MODEL 

CROSS REFERENCE TO RELATED APPLICATIONS 

This application claims the benefit of U.S. Provisional Application 
Number 60/173,602, filed December 29, 1999, which is hereby incorporated by 
reference in its entirety. 

BACKGROUND OF THE INVENTION 

This invention relates generally to predictive modeling, and more 
5 particularly to predicting train consist reactions to specific stimuli. 

Typically, when a train consist moves from one location to another, 
knowledge of the effects of incurred stimuli of the consist as the consist moves along 
the tracks is useful for increasing the efficiency consist movement. Therefore, a train 
consist predictive model would aid railroad transportation companies in improving 
10 operations. However, modeling the movement of a train consist is complex and time 
consuming. 

Therefore, it would be desirable to provide an effective method of 
modeling movement of a train consist so that kinetic characteristics, such as stopping 
distance, speed, and acceleration can be determined and utilized to improve the 
1 5 movement of freight from one location to another. 

BRIEF SUMMARY OF THE INVENTION 

In an exemplary embodiment, an adaptive train model (ATM) includes 
a semi-empirical mathematical model of a moving train consist used for simulating 
and predicting train reaction to external stimuli. The consist includes at least one 
locomotive and at least one railcar. The ATM utilizes a system including at least one 
20 sensor located on the train consist, a database, and a computer. 

More particularly, the ATM is used for predicting such things as, train 
acceleration, train speed after one minute, and a shortest braking distance of the train 
consist. The model is adaptive because it is built and updated while the train is 
moving. The model has a set of unknown parameters that are updated as the train 
25 consist is moving and new data are collected. Then, based on the most recently 
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measured data, those parameters are used for predicting movement characteristics of 
the train consist. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 is a schematic of a system used to implement an adaptive train 
model (ATM) in accordance with one embodiment of the present invention; 

5 Figure 2 is a schematic diagram of a compressed air system, with an 

auxiliary reservoir, used in a train consist utilizing the system shown in Figure 1 ; 

Figure 3 is a chart depicting the difference between actual stopping 
distance and the predicted stopping distance determined utilizing the ATM with no 
initial values; 

10 Figure 4 is a chart depicting the stopping distances for a consist 

determined by the ATM using initial values previously collected; 

Figure 5 is a chart depicting the stopping distances for a consist 
determined by the ATM using information collected for the consist on a prior trip 
where the load was different from the present load; 

15 Figure 6 is a chart depicting variations of the predicted stopping 

distance for a consist determined by the ATM using as initialization information only 
some of the data collected based on values from a prior consist; 

Figure 7 is a chart depicting the predicted stopping distance for a 
consist determined by the ATM using modified initial data; 

20 Figure 8 is a chart depicting the predicted stopping distance determined 

by the ATM for the stopping distance for a consist as the brakes are fully applied 
instantly, using no initialization values; 

Figure 9 is a chart depicting variations of the predicted stopping 
distance determined with the ATM using initialization information based on data 
25 collected after three trips on the same track, and using the predictions made on the 
third trip as references; and 

Figure 10 is a chart depicting an error of an estimated safe stopping 
distance for a consist using the ATM initialized with weighted information. 
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DETAILED DESCRIPTION OF THE INVENTION 

As used herein, the terms "train consist" and "consist" mean a train 
having at least one locomotive and at least one railcar. However, a "consist" may 
include a plurality of locomotives physically connected together with one locomotive 
designated as a lead locomotive and other locomotives designated as trailing 
5 locomotives, and a combination of railcars (freight, passenger, bulk) physically 
connected to the locomotives. 

Figure 1 is a schematic of a system 10 used to implement an adaptive 
train model in accordance with one embodiment of the present invention. System 1 0 
includes at least one sensor 12 located on a railcar (not shown) of a train consist (not 

10 shown), and a computer 14. Computer 14 includes a database 16, a processor 18 
capable of executing all functions of computer 14 and a memory 20 sufficient to hold 
data and algorithms used to implement the adaptive train model. Although system 10 
is shown having only one sensor 12, system 10 may include a plurality of sensors 12. 
Sensor 12 monitors known train consist stimuli relating to moving train consist 

15 kinetics, such as, speed, pressure in consist brake pipe, consist position on a track, 
track configuration, rolling forces, aerodynamic forces, elevation caused forces, 
braking forces caused by direction changes, and braking forces caused by application 
of railcar brake shoes (not shown). Sensor 12 communicates stimuli data to computer 
14, which uses processor 18 to interpret the stimuli data. Processor 18 utilizes a 

20 modeling algorithm, described below and hereinafter referred to as an adaptive train 
model (ATM), to process the stimuli data, thereby producing coefficients and 
operational data that are used simulate train consist kinetic characteristics and predict 
train consist movement. In one embodiment, the coefficients and operational data are 
stored in memory 20. In an alternative embodiment, the coefficients and operational 

25 data are stored in database 16 according the ATM. 

The ATM is a semi-empirical model using a set of initially unknown 
coefficients. Those coefficients are determined and adjusted based on the most recent 
collected stimuli data and train reaction to those stimuli, such as acceleration. Once 
determined, the coefficients are used for predicting train kinetic characteristics for a 
30 period of time. The ATM is semi-empirical because it takes in account the most 
important elements of consist kinetics but considers as unknown their interaction. 
The ATM makes the assumption that the train consist length is constant, the load the 
consist is carrying is uniformly distributed along the consist, the brakes on all railcars 
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have the same unknown efficiency, and the spiral curvature (1/R) of the track has a 
linear variation with distance. 

Table 1 , below, defines the notations that will be used to describe the 
ATM and how the notations are used to predict consist kinetic characteristics. 
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TABLE 1 



Symbol 


Meaning 


Unit 


a 


Acceleration 


feet/s^ 


a(t) 


Acceleration as a lunction ot time 


feet/s^ 




brake tunctions 


feet/s^ 


C,(t) 


braking effect caused by lateral triction when tram is in 
curve 


feet/s^ 




braking el'tect caused by weight increase when train is m 
curve 


feet/s^ 


u 


distance 


teet 




distance as function ot time 


feet 




degree ol a curve (angle tor lUO teet ot track)' 


degrees 


H(t) 


Elevation tunction 


feet/s^ 




horce 


pounds 


g 


Gravitational acceleration (9.82 m/s" = 32.218 teet/s ) 


feet/s^ 




Corrective factor lor the ettect ot aerodynamic triction 


no unit 


Kb, 


brake lunction coetticient 


no unit 




Corrective tactor tor the ettect ot dynamic brake 
application 


no unit 




Corrective factor tor the ettect ot elevation change on 
segment i of the train 


no unit 




Corrective tactor for the ettect ot lateral triction when 
train is in curve 


no unit 




Corrective tactor tor weight increase when tram is m 
curve 


no unit 




Corrective tactor tor friction ot a tram rolling on 
straight horizontal track 


no unit 


K„ 


release tiinction coetticient 


no unit 




Dynamic corrective factor tor triction ot a tram rolling 
on straight horizontal track 


no unit 




Corrective factor tor the ettect ot throttle application 


no unit 


L 


total tram length 


teet 




length ol segment / 


teet 


i„ 


length ot the segment i section J ot the tram 


teet 


M 


total tram mass 


lbs 




mass ot segment / 




m,, 


mass ot the segment i section j ot the train 


lbs 




Number ot powered axles 




p(t) 


Pressure in brake pipe measured at front locomotive 


psi 




Maximum pressure in brake pipe 


psi 


R 


curve radius 


teet 




release tunctions 


teet/s^ 


L 


tram length 


feet 


■i'(t) 


traction force 


pounds 


V 


speed 


feet/s 


v(t) 


speed as tunction ol time 


feet/s 


vd 


speed recorded in database 


feet/s 


W 


total tram weight 


lbs 


w„ 


weight ol the segment / section J ot the tram 


lbs 



1 The field CURVE in track database. 
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Model Description 

The ATM determines the train kinetic elements based on a force 
balance. As the train moves, there are forces pulling the train forward and other 
5 forces pulling the backward. The forward forces accelerate the train movement, while 
the backward forces have a braking effect and create a deceleration. The effect of the 
vertical forces are considered only as an element affecting the rolling friction of the 
consist. 

In an exemplary embodiment, a balance of forces applied to a consist 

10 at time t is 

F(t) = M (K, + K„ v(t))+ K, v(t)' + 
M K„ Ei(t) + M EjCt) + M K,3 E,(t) + M K,, E,(t) + 
MKp Cp(t) + MK, C,(t) + 
Kb, B,(t) + Ki,2 B^O) + 3,(0 + B^(t) + 
1 5 K,, R,(t) + K,2 R^Ct) + K,3 K,(t) + K,4 R4(t) + Kj D(t) + K, T(t). (I) 

The meaning of the function and coefficients is described below. 

A rolling forces term, M (K, + K„ v(t)), is the sum of all forces that act 
on a train moving along an ideal straight horizontal line when no brake or throttle is 

20 applied. Those forces can be determined on each independent railcar. They are caused 
by the friction between rail and wheels, rail elastic deformation under railcar weight, 
deformation (compression) of the rail understructure caused by railcar weight, bearing 
friction, and other forces in direct relation with the railcar weight and railcar speed. 
Adding the effect of all railcars, gives a force that is dependent on train total mass and 

25 speed. 

K, is the only coefficient in the ATM that does not have an associated 
element. K^^ is associated with train speed. As the consist speed changes, K„v 
modifies the value of the rolling forces. 

Aerodynamic forces are expressed by the term, v(t)^. The 
30 aerodynamic force's are caused by the friction between the moving train consist and 
the stationary air. The total force depends on locomotive aerodynamic profile, railcar 
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aerodynamic profiles and the sequence of different railcar types along the consist, 
such as, flat bed, single stacked, double stacked, caboose, gondola, etc. For such 
reason the value of the coefficient should be different from one consist to another as 
each consist has its own aerodynamic profile. 

5 Elevation caused forces are expressed by M (K^, Ei(t)+ K^2 ^li^)^ 

£3(1)+ E4(t)). It is assumed that the consist has up to four segments, each 
uniformly loaded. This assumption is important for modeling the movement of very 
long consists, where distribution of weight has an important contribution on the speed 
determination. The size of each segment is determined from loading information 

1 0 provided in a consist description. If such information does not exist, all four segments 
are considered to be of equal length. For small consists (less than 20 railcars) or 
uniformly loaded consists, such as coal, ore or other specialized consists, the whole 
train is considered as one segment. 

As a train travels over hills and valleys, some sections of the train are 
1 5 moving uphill, while other train sections are rolling downhill. 

A force f,j created by section j of segment i of the consist is tangential 
to the ground and has the value 

fy = w,j cos(7X /2-a,j). (2) 

20 Factor w^ is the weight of consist section] of segment i. Since the load 

is assumed to be distributed uniformly along the segment, the weight of the section j 
of segment i should be 

w,j = M^gl,/Lj (3) 

25 Then an elevation effect on consist section j of segment i is 

fy = M, g l,j / L, cos(7i /2-ay). (4) 

Values l,j are the "footage" determined by a track database using two 
successive elevation records. The first and the last sections of the consist should be 
30 determined according to the position of the head and the tail of the train with respect 
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with the surrounding terrain. As the train travels, the size of the train sections are 
changing continuously, thus they are time dependent. Constraint 

is verified at each moment t. The value of cos( 71 / 2 - cQ is approximately the same 
as the value of the field "slope" in a track database divided by 100 (PLEASE 
DEFINE FIELD SLOPE). 

E,(t) is determined by 

E,(t) = l,j(t) slope, / 100 DB_DIRECTION. (6, 

The DB DIRECTION factor defines the vertical direction each 
segment of the consist is moving. If the segment is on a downhill slope, 
DB_DlRECTION is set to negative one (-1), otherwise its value is one (1). 

Term g M, /L; is a constant that was not included 'm(6). Thus, value of 
Kg, should be g M, /L„ giving 

K,. = gM,/L, (7, 

From equation (7) 

K,,L, /g = Mi^ and fS 
4 

g^K„L, = Mg = W. (9 

Noise associated with the data collection and database inconsistencies 
are not considered, therefore, the sum is an approximation of train weight. 

Braking forces caused by direction change are expressed by the term 
M (Kp Cp(t) + K| Cj(t)). Each change in track direction has three distinctive segments, 
the entry easement, the circular curve and the exit easement. The entry easement is a 
spiral with the curvature varying from 0 to a final value of 1/R. In other words, the 
radius of the turn varies along the easement from 00 to R. The curve is a circle with a 
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constant radius of R. The exit easement is symmetrical to the entry easement. It starts 
with radius R and ends with a radius equal with oo. When a railcar is moving at 
constant speed along the circular section of the track, it is affected by the centrifugal 
acceleration determined by 

F,= mvVR. (10) 

When the railcar is moving along a straight track, 1/R is zero and the 
centrifugal force is nonexistent. The reason for building the entry end exit easements 
is to prevent the centrifugal force discontinuities from zero to the value created by the 
movement on circular track. Radius R is determined by using a field from track 
database, giving 

R = 180/71 100/Dc= 5729.58 / Dc (11) 

The centrifugal force is compensated for by raising the outer rail of the 
track, thereby tilting the railcar inward. The point of the railcar center of gravity O 
depends on the level of railcar loading and the type of load. There are two main 
forces acting on a railcar moving on a circular track, a railcar weight (w) and a 
centrifugal force (F^,). The centrifugal force is dependent on railcar mass, speed and 
track curvature. The railcar weight is dependent on railcar mass and gravitational 
acceleration. These two forces are composed into a resultant (w'). In order to 
prevent train derailment, the position of this force should always be between the two 
rails. Therefore, 

w' = sqrt{ + F/), and (!2) 
w' = m sqrt(g2 +ivVRf) = m sqrt( (R g)'+v')/R. (13) 

The ideal position of this force is perpendicular to the track. Such 
position corresponds to a speed vd recorded in the data base for that particular curve 
and is strictly related with the amount of superelevation. The force w' is decomposed 
in two components, one is perpendicular to the track (Fp), while the other is parallel to 
the track (F,). The component perpendicular to the track is the "apparent" weight of 
the railcar and the component parallel to the track creates lateral friction. Therefore, 

Fp = w' cos(e) = m sqrt( (R g)'+v") cos(9)/R, and (14) 
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F, = w' sin(e) = m sqrt( (R g)'+v'') sin(e)/R. • (]5) 

The angle 0' is between forces w and w', giving 

tan(9') = F,/ w = ( m vV R)/(m g) = vV(R g). (16) 

@ is the superelevation angle and should be equal to the angle between 
forces Fp and w. The superelevation angle should satisfy the relation 

tan(0) = vdV(R g), (I7) 

sin(©) = vd" /sqrt( (R gf + vd"); cos(0) = R g/sqrt( (R g)' + vd"), and (18) 

sin(e) = /sqrt( (R g)^ + v^); cos(e) = R g/sqrt( (R g)" + v^. (19) 

The angle 0 is the difference 0' - ©. Therefore, 

cos(e) = cos (9' - 0) = cos(e') cos(0) + sin(e') sin(0), and (20) 

sin(e) = sin (9' - ©) = sin(e') cos(©) - cos(0') sin(0). (21) 
Which gives 

cos(e) = ((R gf W vd^) / ( sqrt{ (R g)^ + vd^) sqrt( (R gf + v^)), and (22) 

sin(e) = R g (v^ - vd^) / ( sqrt( (R g)^ + vd^) sqrt( (R g)^ + v^)). (23) 
Replacing the values in (14) and (15), gives 

Fp = m ((R g)" W vd') / (R sqrt( (R g)' + vd^)), and (24) 

F, = m g (v^ - vd^) / sqrt( (R gf + vd*). (25) 

These two forces are associated to the sections of track that are not 
tangent (straight). Fp represents the braking force caused by the effect of gravitation 
and centrifugal acceleration. If the effect of gravitation was already included in 
rolling forces, the relation in (24) should be adjusted by subtracting m g, giving 

Fp = m (((R gf W vd^) / (R sqrt( (R gf + vd^)) - g). (26) 
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The effect of all sections of the train that are on a circular track at a 

moment t is 

Fp = Z m, (((R, g)^ vd I (R, sqrt( (R, g)^ + vd - g), (2 7) 

5 and using (15) gives 

F, = |j m, g (v^ - vd / sqrt( (R, g)^ + vd Z). (28) 

Value m, is the mass of the train on each circular section i of track. For 
a uniform train loading the value should be proportional with the section length. The 
10 final expression of forces is then given by 

Fp = M/L Z 1, (((R, g)^ W vd I) I (R; ^ sqrt( (R, g)^ + vd - g), and (29) 

F, = M/L S 1, g (v^ - vd ,2) / sqrt( (R, g)" + vd (30) 

Similar formulae are associated with forces developed by the train 
1 5 sections in easements. Since the radius of easement is variable, forces associated with 
easements are considered to be half of the forces created in circular sections. Such a 
supposition is correct if easements have a linear variation of a curvature 1/R. 
Removing constant terms, the expression of the "perpendicular force effecf for 
sections of train located in curve sections of track is 

20 Cp(t) = Z l,(t) 5, (((R, gf +v^ vd / (R, sqrt( (R, gf + vd ,^)) - g), (3 1) 

and the "lateral force effecf is 

C,(t) = Z l,(t) 5, (v^ - vd / sqrt( (R, g)^ + vd ,^)) . (32) 

25 where 5; are 1 for circular sections and .5 for easements. 

Braking forces caused by the operation of the consist brakes are 
expressed by B,(t)+ B2(t)+ Ki,3 B3(t)+ B4(t). The sum of four terms K^; B,(t) 
denote the total braking force applied to the train. The value of this sum will differ 
-11- 
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from zero only when the brake shoes are pressing the wheels on at least one railcar. 
After an engineer reduces the pressure in brake pipe, the brake shoes along the consist 
are successively applied against the wheels increasing the total braking force to a 
maximum value, railcar by railcar. Although all shoes are pressing the wheels with 
5 the maximum force, the total braking force does not remain constant. As the train 
speed decreases, the friction coefficient between the wheels and shoes changes its 
value modifying the value of the frain total braking force. 

The friction coefficient between the wheels and shoes is in direct 
relation to the material composition of the brake shoes. Brake shoes can be made 
10 from cast iron, cast iron with high phosphorous, or high friction composition blend of 
organic and inorganic materials. In addition, the static and dynamic friction 
coefficients are different from one shoe type to another. The following formulae 
determine the friction coefficient. 

For cast iron shoes the friction coefficient is, 
15 C, = .18+.28*exp(-.25*v), and (33) 

For high composition shoes the friction coefficient is, 



C, = .255+.ll*exp(-.07*v). (34) 

20 Figure 2 is a schematic diagram of a compressed air system 30 for a 

train consist using an auxiliary reservoir. Compressed air system 30 includes a 
consist brake pipe 32 that runs the length of the train consist, a main reservoir 34 for 
storing pressurized air, an auxiliary reservoir 36 for storing auxiliary pressurized air, a 
brake cylinder 38 that engages and disengages the brake shoes of a railcar (not 

25 shown), and an exhaust 40 that allows air to be vented from system 30. Although 
system 30 is shown having only one brake cylinder 38, each wheel of a railcar is 
coupled with a brake cylinder 38 such that a respective set of brake shoes are operated 
for each wheel. Brake application is caused by a reduction of air pressure in consist 
brake pipe 32. The dragging force created by the brake application is equal to the 

30 product between a friction coefficient C/y) and a force f^,(^) applied by the brake shoes 
to the wheels of each railcar, as the pressure builds up in brake cylinders 38. The 
value of the pressure in brake cylinders 38 relates to a pressure drop in brake pipe 32 
Ap and an initial pressure in the auxiliary reservoir 36 Pa. 
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The variation of pressure in brake cylinders 38 depends on a variety of 
factors, such as, the length of brake pipe 32, air temperature, humidity, the operational 
status of a brake valve 42, and the amount of slack in the brake shoes. However, there 
are likely more factors, some unmeasurable or unknown, that effect the pressure in 
5 brake cylinder 38. To compensate for variations two models that evaluate the brake 
cylinder force in two opposite conditions are used, one supposes that the pressure is 
building fast in brake cylinders 38, while the other considers that the pressure is 
building much slower. Although, the two pressure models are designed to match the 
braking force developed by short trains and long trains, their linear combination can 
10 mimic the time variation of the braking force for a large variety of trains in many 
external conditions. The proper combination of those two pressure models with the 
friction coefficients is a critical task of the Adaptive Train Model. 

The model for fast building pressure uses a pressure drop causing 

equalization 

15 Ep = (Pa "^285856 + g4_|49g89g9)/ 12.43010076, (35) 

a reference pressure 

Rp = min(Ep, Ap), and 

a rigging equation 

20 Rigg = min(0, max(l, (12.4301008 * Rp - 84.1496897) / Pa' 's^s'sse^). (36) 

The rigging value should be between 0 and 1. Therefore, brake 
cylinder maximum force equals 

BCf = (345.0445077 * Rp - 4986.482562) * Rigg, and (37) 

25 

modified time equals 

T = max(0, (t -2.5332561 1) * (297.383526 / (0.08709894 * L +1)) <'™'^'<^ 
* [ 1 / (-0.02072556* Pa- 0.00036791 * Pa' + 6.47422981) 
+ 1 /(-0.00423 106 * Rp + 0.00 137952 * Rp' +0.02506324)]. (38) 

30 

If Ap > Ep, the modified time is replaced by 
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J = J (1+25 2474268 (Rp - Ep)/Pa -1 48 737743 ((Rp - Ep)/Pa)2 ) ^g-^ 

Therefore, the braking force for fast building pressure is given by 

Bff = min(0, max(l, (T + 3.86950758 * V + 0.23164628 * T') / 
(16367.9101 + 1 1 1.652789 * T + 27.6134504 8 * - 0.0026229 * T^) )) BCf (40) 

The model for slow building pressure uses a pressure drop causing 

equalization 

Ep = (Pa ' •s«'22«6+ 83.5394856) / 12.47984, (4!) 

a reference pressure 

Rp = min(Ep, Ap), and (42) 

a rigging equation 

Rigg = min(0, max(l, (12.47984 * Rp - 83.5394856) / P, ' (42) 

The rigging value should be between 0 and 1. Therefore, brake 
cylinder force is given by 

BC3 = (346.46923 * Rp - 5015.0018) * Rigg, and (44) 

modified time equals 

T = max(0, (t - 4.99998034) * [297.383526/ (0.08709894* L +1)] " 
* [ 1 / (-0.00781443 * Pa - 0.00034658 * Pa" + 6.76259649) 

+ 1 / (0.07532957* Rp + 0.0020505 * Rp' + 0.93285061) ]. (45) 

If Ap > Ep, the modified time is replaced by 

7 ^ = (1+2 47330639 (Rp - Ep)/Pa - S 9273390S ((Rp - Ep)/Pa)2 ) ^^^^ 

Therefore, the braking force for slow building pressure is given by 

Bf, = min(0, max(l, (T, + 2.00986206 * T,' + 0.81412194 * T/) / 

(0.00067603+ 169.361303* T, + 8.95254599* T,^ + 0.58477705* T^^) )) Be, (47) 
-14- 
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Taking into consideration the model for fast building pressure and the 
model for slow building pressure, the total braking force has four components 

B, = Q Bff , 
B2 = C,Bff, 
B3 = C, Bf, , and 
B4 = C,Bf, 

The balance of forces equation (1) also takes into consideration the 
brake release pressures. When an engineer restores the pressure in the brake pipe to 
the initial value, the brake shoes are removed from wheels and the dragging force 
created by the brake application decreases. The ATM considers brake release as a 
force opposing the braking force. The ATM has to identify the proper coefficients 
such that the release force cancels the brake force when shoes are in fact completely 
removed from wheels. As in the brake application, two conditions are considered, a 
fast release model and a slow release model. 

The model for fast release utilizes a pressure drop causing equalization 

Ep = (Pa • + 84.14968969) / 12.43010076, (48) 

a reference pressure 

Rp = mm(Ep, Ap), and (49) 

a rigging equation 

Rigg = min(0, max(l, (12.4301008 * Rp - 84.1496897) / '^^^sse^-j^ (^q) 

The rigging value should be between 0 and 1. Therefore, brake 
cylinder maximum reference force is given by 

BCf = (345.0445077 * Rp - 4986.482562) * Rigg. (51) 
Thus, the fast release force is given by 

Rff = min(0, max(l, (t + 3.86950758 * f + 0.23164628 * t') I 
-15- 
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(16367.9101 + 111. 652789 *t + 27.6134504 8 *t'- 0.0026229 *t')))BCf (52) 
The model for slow release utilizes a pressure drop causing 

equalization: 

Ep = (Pa 1 83.5394856) / 12.47984, (53) 

a reference pressure: 

Rp = min(Ep, Ap), and 

a rigging equation: 

Rigg = min(0, max(l, (12.47984 * Rp - 83.5394856) / P, ' '^^^zm^). (54) 

The rigging value should be between 0 and 1 . Therefore brake cylinder 
reference force: 

Be, = (346.46923 * Rp - 5015.0018) * Rigg. (55) 
Thus the slow release force is given by: 

Rf, = min(0, max(l, (t + 2.00986206 * + 0.81412194 * t') / 
(0.00067603+ 169.361303* t + 8.95254599* + 0.58477705 * t') )) Be,. (56) 

Taking into consideration the model for fast release and the model for 
slow release, the total release force has four components 

R, = C,Rff, 
R2 = C,Rff, 
R3 = C, Rf, , and 
R4 = C,Rf, 

The dynamic brake force is represented by D(t). The value in 
pounds of D(t) is determined by 

D(t) = Fv"N3,/591.43, (57) 
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where I is a current in Amperes, v is a speed in feet/s and is a number of powered 
axles. 

The traction force is represented by K^ T(t), where T(t) is the traction 
effort and is dependent on locomotive type and train speed. K, will be one if the train 
has only one locomotive, the locomotive efficiency is 100% and the locomotive 
wheels do not slip. 

For trains using more than one locomotive type, or not using the same 
number of locomotives all the time, T(t) will be the sum of traction efforts developed 
by ail active locomotives. 

Since a train mass M does not change value during the trip, the relation 
in (]) is divided by the total train mass. The resuh is the expression of train 
acceleration 

a(t)=K, + K„ V + K, + K„ Ei(t) + K,^ E^Ct) + K^s E,(t) + K,, E,(t) + Kp Cp(t) + K, 
C,(t) + Kb, B,(t) + Kt,, B^Ct) + B,(t) + K,4 B^Ct) + K„ R,(t) + K,, R,(t) + (58) 
K,3 R3(t) + K,4 R4(t) + Kd D(t) + K, T(t). 

In this relation all current and past values of speed, acceleration, 
throttle position, pressure in brake pipe, and position on track can be used for 
determining coefficients K,. Once the coefficients are determined, they can be used 
for predicting the consist speed, acceleration, and other values at any moment of time. 

Determining the Coefficients. 

In one embodiment, to determine the coefficients for the ATM, the 
Least Squares Method (LSM) is utilized. In one embodiment, the behavior of a real 
world object depends on a set of measurable independent stimuli x,, Xj, . . . x„. To any 
set of stimuli, the object provides a measurable reaction y. Therefore a mathematical 
relation that approximates the reaction of the object to those stimuli is developed. 

The independent set of variables x (stimuli) and the dependent variable 
y (reaction), are measured. Measuring those variables, m times gives the following 
table of measurements M: 
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X2,i ^2,2 ■ ■ ■ ^2,n 

X^,, ■•• ^m,n Ym- 

The dependent variable y can be determined by a linear combination of 
other independent variables x, giving the equation 

y = ZajXj. (^9) 

The dimension m of the table is variable. At any time another 
experiment can be made and the new set of values added to a previously buih table. 

For any set of chosen or determined values of coefficients an 
approximation y' is determined. The difference between the measured y and the 
10 determined y' is the approximation error s, where 



= y.- ?i 



X for i = 1, 2, n. 



By choosing different sets of a^ coefficients, the errors b, change. In 
one embodiment, the values of aj that minimize the expression 



E= S 8,^= Z(y,-E a^x,J^ (61 

are calculated. The global error E depends on values only, while all other elements 
in relation (61) are measured numerical values. Then the partial derivative of E with 
any a,, will be zero, for example, 

SE/5ak=0fork=l,2, ... n. (62, 

Therefore, replacing the value of E from (61) in (62) and 
differentiating gives the equation 

-2 £ ((y, - aj x,_ J ) (d aj /8 a, x,_ j ))= 0 for k= 1 , 2, . . . n. (63 

The term d Id a^ is 1 when j is equal with k and zero in all other 
cases. Ignoring the constant term -2 gives the equation 
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S((y,-S ajX,,j)x,,,) = 0fork=l,2, ... n. (64) 

Changing the summation order and moving known elements on the 
right hand of the equation gives 



Z aj(E x,,i x,J=E y,x,,, fork=l,2, ... n. 



(65) 



Both summations by i are using measured numbers. Those 
summations provide the matrix that is used for determining the unknown a, 
coefficients. Using the notations 

=Xj k= 2 x, j X, k, and 

Yk= E y, x,,k, (66) 

gives the final form of the system of equations 

Z Xk,a,=Ykfork=l,2, ... n. (67) 
j=l 'J 

The relation (67) is a system of n equations with n unknown variables 
a,. The matrix of this system X and the right hand term Y are determined from the 
measured data. The values of unknown coefficients aj are found by solving the system 
of equations (67). After determining coefficients a.^, the individual error for each 
measurement is determined using 

s, = y, - S a, X, J for i=l , 2, . . . m. (68) 

Replacing the expression of error C; in (65) gives 

2s,x, k=0fork=l,2, ... n. (69) 
If one of the x, j has the value one, then the individual errors also have 

the property 

S e, = 0. (70) 
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If measurements are not all made in the same conditions, weighting 
factors reflecting our confidence in the validity of each measurement are used. In 
such a case, the expression is minimized by 

E = E w. s,^ = S w, (y - aj x,_ j )^ (71) 
where each independent error has a different weight. Additionally, 

Xk,j =Xj k= S w, X, J X, k, and 

Y,= Z w,y,x,,,. (72) 
If the approximation has a free term, the property 



S w, s, = 0 appHes. (73) 

The relation (58) defines the train acceleration. The following table is 
a mapping of the train movement model to least square method notation. 



i'rain 






Ao 


1 


Xo 




A: 


V 


X, 





A, 








A3 


K,(t) 


X3 








X4 




as 




X5 




a6 


Jti4(t) 


Xg 




% 




X7 




a^ 


C,(t) 


Xg 


Kb. 


a? 




Xg 
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Train 


im 


^b2 




B2(t) 


X,o 


^b3 


Bij 


B3(t) 


^11 




^12 




X]2 




a,3 


RiCt) 


X|3 






R^Ct) 






ai5 


RsCt) 


X,5 






R4(t) 






ai4 




Xi4 




aj7 


■i(t) 


Xj7 


acc(t) 





The measured parameter Xg always has the value 1. Time t is not 
explicitly used in this model. All x, values are determined using the latest measured 
acceleration, speed, train position and pressure in equalizing reservoir 36 shown in 
5 Figure 2. 

The X matrix is symmetrical with the positive diagonal. As a result, 
indefinitely adding elements to the matrix causes data degeneration and loss of 
precision. In order to prevent data degeneration each measurement is weighted. 
Thus, the weight of the present measurement is set to one, the weight of the 
1 0 measurement before that is set to a smaller value and the weight of each preceding 
measurement is set to yet a smaller value. 

The value s is set to a positive number smaller than 1 and the series 1 , 

(1-s), (1-s)^, (1-s)^ . . . converges to zero. By considering the series of weights for 
consecutive measurements, after an infinite number of measurements, any positive 
15 term on the diagonal of the LSM matrix will have a value limited according to the 
equation, 

x[i,i]<3^|^ (l-s)''=3^(l-(l-sr)/(l-(l-s)) = 3Vs, (74) 
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where 3 is the maximum value between one and the largest absolute value of the 
considered parameter. From (74), applying the weighting procedure presented above 
to LSM generates a matrix that has all terms on the diagonal smaller than 3^ /s, where 
0 < e < 1. Additionally, equally weighted measurements (corresponding to s = 0) 
create a condition for floating point errors, and that large values of 8 (almost one) 
assure better precision for number representation, but gives too small of a weight to 
prior samples, practically eliminating them from selection. 

The optimal value will be somewhere between zero and one. Using s = 
10 with 32 bit floating point representation provides 4 correct decimal digits. Using 
s =10'^ with 64 bit floating point representation provides 1 1 correct decimal digits of 
the result. 

The term x[l,l] is associated to the free term. Without any weighting, 
this number is equal with the number of samples considered for determination. When 
using weights, the term is called "the apparent number of samples". Using the 
weighting procedure the value of x[l, 1] converges to 1/e. 

The sum in (74) is performed for an infinite number of terms and there 
is no need to store more data than the current evaluated x vector. The weighting of 
data is realized by "aging" the X matrix and the Y vector before adding the latest 
collected data. If X(n) is the matrix X for step n and Y(n) is the vector Y in the same 
step, the next step elements are 

X(n+l) = (l-8)X(n) + ||x„„^j x„„,,||forj,k=l,2,3, ...m , and 

Y(n+1)= (1-s) Y(n) + y„+| ii x„+, , || for k=l, 2, 3, ...m. (75) 

From (75), the following corollary is derived. 

IfX(n) has the inverse P(n) andX(n+l) =X(n) +x(n+l) x(n+l) T then 
P(n+1) = P(n) - [P(n) x(n+l) ] [[x(n+l)T P(n) J / (x(n+l)T P(n) x(n+l) + 1)]. (76) 

However, if there is problem in the initial matrix, X(0) is null and it 
has no inverse. To overcome this inconvenience without destroying the precision the 



-22- 



GEH-0 1-064 



initial value of X(0) is set to some small values on the diagonal, for example 10"', with 
the rest of the matrix set to zero. During the process the small value is overwritten by 
larger x(n)x(n)^ values. The inverse of X(0) is P(0) with all elements set to zero 
except for the diagonal that has the inverse of each element from the diagonal of X(0) 
5 divided by m. 

A Y vector is updated according to 

Y(n+ 1 ) = Y(n) + acc x(n+ 1 ), and (77) 

the values of aj coefficients are determined according to 

10 a(n) = P(n) Y(n). OS) 

When performing the matrix operations the total number of computer 
operations is 13m^+2m+l. The Gauss- Jordan triangulation method, requires for 
solving the same problem, m^+(m^-m)/2 operations. 

15 Model Confidence 

All variables in the ATM are independent. Therefore, all coefficients 
determined will be independent. Thus, the real train behavior and measurements are 
not affected by large errors, and after a short period of time, the value of each 
coefficient will be approximately the same from one second to another. 

20 The stability of coefficient values effects the stability of predictions 

that are provided. One way to measure the quality of information gathered by the 
model is to monitor the differences between the predicted speed and the measured 
speed. The speed may be predicted for one minute, 15 seconds, or only one second. 
The potential disadvantage is that the quality of a prediction at a time after the 

25 prediction was actually made is unknown. 

For a one second speed prediction the response is available after one 
second, but is affected by the measurement noise, and the one minute prediction 
provides a response after one minute, which is too late. A five second or fifteen 
second prediction is sufficient for judging the quality of model predictions. 
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Each speed measurement is associated with the value v^, predicted a 
few seconds earlier. For a large set of measurements or predictions the variance a is 
associated according to the equation, 

<j^= |j(v,-v,)V(n-l). (79) 

The only values that have to be stored are n, the sum of the square 
values ( Vj,)^ and the sum of the values v^. As a new value v^ is predicted, it is 

stored in a circular array from where it is retrieved when the corresponding V;, 
measurement is available. The differences between the values are used for updating 
the value of the variance a^. The confidence coefficient is determined as a function 
of \(vk-v]^/o\ or fv^v/t//a^in accordance with the following table. 
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48.4 


1.0 


1.0 


31.74 
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2.3 


13.36 


2.0 


4.0 


4.56 


2.6 


6.J 


1.24 


3.0 


9.0 


0.26 



The confidence column is determined by the modified ERF (PLEASE 
DEFINE ERF) function, 

; 1 J _l2 

V 27j: zs_ 

a 

Value s is the absolute value of the difference between the last 
measured value and the corresponding predicted value. The relation. 
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1/(1+0.77439865 x + 0.80155580x2 + 0.54110010x3 + 0.030646485x4), 

where x is le^/cr |, provides an approximation of the ERF function with errors between 
-.00067 and +.001315 for 0 < x < 3.4. 
Predictions 

The ATM determines the instantaneous train consist acceleration. 
Predicted values such as acceleration, speed, and distance are integrated starting with 
the present position, speed, and acceleration. If all the coefficients are already 
determined, the equation (58) is written as 

a(t) = 0(t). (80) 

The beginning of the prediction is considered as t=0. The numerical 
integration is performed in steps. Thus, step i=0 corresponds to t=0, step i=l 
corresponds to t-At, and so on. The size of At is capable of being equal to one 
second, a half second, .1 seconds, or smaller and larger values. A better prediction 
will be provided using a smaller value so that the value ao equals the acceleration 
value measured at t=0. 

When determining the acceleration, the <l>(t) function uses two 
determined values, the current locomotive position and the speed. To maintain 
consistency with other relations, the function parameter used is i. The notation 0(i) 
represents the value of the function determining the acceleration, using a speed value 
Vj and a distance dj. 

In one embodiment, the speed is determined by utilizing the known 
Euler method, giving the equation 

v,=v,.,+ Ata,_,. (81) 

The distance is determined by the speed value at the beginning and the 
end of the time interval, according to 

d, = d,.i + At (v, + v,.,)/2. (S2) 
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The acceleration is provided by the relation (58) symbolized as the 
0(i) function. 0(i) uses the new position d, to identify the length of the consist 
affected by elevation change and track curvature and the most recent determined V; 
speed, accordingly 

a, = 03(1). 

The acceleration is considered to be constant for the length At. The 
error in determining the distance d„is approximately At ' 0(t) d,/12 where 0 < x < nAt. 
When (D(t) has negative values (braking) a shorter distance is given. Therefore, using 
the Euler method for approximating the stopping distance, provides smaller values 
than actual stopping distance. 

The relative error of distance determined with the Euler method is 
proportional with At\ and requires one evaluation of function 0(i) per iteration. 

In another embodiment, the speed is determined by using the known 
Milne method for second degree equations =f(x, y,y'). Each integrated parameter 
is determined twice. A predicted value is determined first and a corrected value is 
determined second. Thus, the speed prediction is given by, 

v,= v,.4+ 4At(2a,.3- a,.2+2a,.,)/3, 
d, = d,.3+3(d,_]-d,.2) + At"(a,_,- a,_2), and (84) 
a,= 0(i). 
The correction is given by, 

v,= v,_2+ At(a,_2+4a,., + a,)/3, 
d, = d,., + At(v,+ v,.,)/2-At'(a, -a,_,)/ 12, and (85) 
a,= a>(i). 

The Milne method requires two determinations of <l>(i) function per 
iteration and the relative error is proportional with At\ 
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Initial values, a^, Vq, do, are the last measured values. Additionally, 
Milne's integration method requires a.j, v.,, d.,, a.j, v.2, d.j, a.3, v.3, d.3, which are the 
values measured in the last three steps before the integration starts. Since the data 
collection time interval is one second, At is one second. 

5 The speed after one minute is determined by using the integration 

procedure presented in (83). During the integration, the function (I)(i) is evaluated 
keeping constant the pressure in the brake pipe and the traction force. Each step of 
integration determines a„ v„ d, from previously determined four triplets. The 
integration continues for sixty steps giving the predicted value of v^q. 

10 The shortest breaking distance is determined by using the same 

integration method until the speed v„ becomes zero. During the integration the 
pressure in brake pipe 32, (shown in Figure 2), is set to a constant value 
corresponding to full service brake application and the traction force is set to zero. 
The integration continues until v„<0, thereby giving the predicted braking distance 

1 5 having a value of d^. 

Model management 

Using information from a consist description, default or initial values 
for model coefficients are determined. The biggest disadvantage of this method is that 
the errors caused by incorrect or incomplete consist information, do not decrease 
20 while new data is collected, so the quality of predictions does not improve with time. 

The ATM uses the system of linear equations given above. The system 
of equations is updated every second with new data about the consist performance. 
Therefore by solving the system of equations an updated model is given every second. 

The normal least squared matrix (LSM) generated for a consist is 
25 suitable to be transported and used as initial data on another consist with identical or 
close characteristics. Information about differences in consist structures is used for 
reshaping the initial normal matrix. 

Figure 3 is a chart, based on an exemplary embodiment, depicting the 
difference between actual stopping distance and the predicted stopping distance 
30 utilizing the ATM with no initial values. The exemplary embodiment relates to a 
consist having two locomotives, 50 loaded railcars and 51 empty railcars, the stopping 
-27- 
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distances are determined by using no initial information. The most important event 
regarding train movement and data collection are given in the following table. All 
determinations refer to the same stopping operation. The X axis presents the time or 
the "age" of the model. The Y axis presents the predicted stopping distance compared 
5 with the actual stopping distance. Each prediction is made using the ATM 
coefficients resulting from data collected prior to the moment of prediction. 



No 


Time 
(sec) 


Trip 
(^leei; 


Event 


Comments 


{ 




0 


"Start nn 

on no c 


y ^.y^ ' FT 1~" ^1 

heavy and moves backward 


2. 


36 


-62 


N otch set to 2 


I'ram keeps moving backward 


3. 


66 


-128 


Notch set to 3 


i'ram starts moving forward after 
backing up 148 feet. 


4. 


126 


-123 


7 psi brake pipe 
reduction 


Tram Speed is about 1 MPH 


5. 


176 


-76 


I rain stops 


irain moved forward tor 71 feet 
of the 148 feet backed up. 


6. 


186 


-76 


K.elease brake 




7. 


186 


-76 


Notch changed from 3 
to 5 


I'ram starts moving forward. 


8. 


685 


8408 


Notch set to 0 




9. 


685 


8408 


7 psi brake pipe 
reduction ^ 


Speed IS 17.i6 MFH 


10. 


716 


9184 


lu psi Dr3,Kc pipG 
reduction 




11. 


746 


9808 






12. 


746 


9808 






13. 


946 


13442 


Notch set to 0 




14. 


946 


13442 


7 psi brake pipe 
reduction 




lb. 


976 


13862 


15 psi brake pipe 
reduction 




16. 


1007 


14076 


i rain stops 




17. 


1176 


14076 


Notch changed from 0 
to 2 


Brakes are still on 


18. 


1176 


14076 


Brake release 


The pressure starts to build in 
brake pipe. 


19. 


1180 


14078 


Train stops 


1 rain moves 3 teet and stops 
because brakes are not released 
yet. 


20. 


1184 


14078 


I rain start moving 


Brakes are releasing gradually 


21. 


1206 


14088 


Notch changed from 2 
to 8 




22. 


2006 


45382 


Notch set to 0 


Speed IS 47.68 Mm 


23. 


2006 


45382 


Brake pipe reduction ot 
26 psi (txill service 
application) 
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No 


lime 


Trip 


Event 


Comments 




(sec) 


(feet) 






24. 


2109 


49609 


Train stops 


Stopping distance is 4227 feet. 



In accordance with the above table. Figure 3 shows that before the 
second 1400, the predicted stopping distance differs from the real value by more than 
+25%. Starting with the second 1447, after traveling 20336 feet, from the previously 
5 collected data, the model predicts the stopping distance that will occur at second 2008, 
with an error smaller than +5 %. After the second 1521, the stopping distance is 
predicted with an error smaller than +2%. Between time 1850 and 2000, the error of 
predicted stopping distance has variations between +5% and +12%. As the predicted 
event become closer, the error of the prediction becomes smaller than +10%. As the 

10 full service brake is applied, the error drops to -4 %. Those variations in prediction 
accuracy are caused by the fact that the ATM determines the best approximation for 
the already collected data corresponding to past events, with the assumption that such 
data are representative for future events. When passed events do not provide enough 
information or are not similar with the future events, the prediction is not very 

15 accurate. 

Analyzing the variation of all model coefficients determined between 
time 1400 and 2000, it is seen that predictions made between seconds 1520 and 1850, 
which are closer to 100%, are made using a set of coefficients that do not properly 
predict the variation of acceleration when the railcar brakes are fully applied at second 
20 2006. Although the acceleration variation is not perfectly modeled, the estimation 
errors compensate for each other and the whole prediction provides a very accurate 
stopping distance. 

Figure 4 is a chart, for a second exemplary embodiment, showing the 
stopping distances determined on a consist having two locomotives, 50 loaded railcars 
25 and 51 empty railcars, using information collected in the previous embodiment. Thus, 
the ATM is not started from zero. The accuracy of the predicted stopping distance for 
the fill! brake application at second 2006, is shown in Figure 4. 

Since the model already has data about stopping the train with a full 
service brake application, the stopping distance is determined having an error between 
30 -4.2% and -2.7%. 
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Thus, the data collected by a train consist during one trip can be used 
as initial data for the ATM when the consist makes another trip with the same amount 
of load. For example, consists moving ore or coal perform the same trip again and 
again with the same consists, loaded with almost the same amount of ore in one 
5 direction and always empty in other direction. 

Figure 5 is a chart, for a third exemplary embodiment, depicting the 
stopping distances determined on a consist having two locomotives, 50 loaded railcars 
and 5 1 empty railcars, using information collected for the consist on a prior trip where 
the load was different from the present load. Figure 5 presents the variation of the 
10 predicted stopping distance when the model is initialized with data collected by 
another train that is 5% heavier and its measured stopping distance is 35% shorter. 

Despite those differences, the stopping distance is predicted with errors 
only between -4% and -10%. This range of errors is caused because although some of 
the initialization information was useful, it also included data that does not fit the 
1 5 characteristics of the consist being characterized. 

Figure 6 is a chart for a fourth exemplary embodiment, based on a 
consist having two locomotives, 50 loaded railcars and 51 empty railcars. Figure 6 
presents the variation of the predicted stopping distance determined with the ATM 
using as initialization information only some of the data collected based on values 
20 from a prior consist. More particularly, only data collected during the 300 second 
period that includes a full service application were used. Thus, the model makes very 
good predictions of the stopping distance after collecting data about the new consist 
for less than 10 minutes. 

When the brakes are applied at second 685, the new information 
25 associated with this event helps to improve the quality of prediction of the stopping 
that is about 1400 seconds away at the time. 

Looking at the variation of the prediction after second 2000, when the 
full service brake is actually applied, the quality of prediction first goes down from 
100% to 96%, then, as the model acquired new data about train reaction to a full 
30 service brake, the quality moved up to 99%. Such variation is characteristic to models 
that improve the quality of coefficients. Since the previous prediction was 100%, it 
means that the global prediction was good, while the acceleration variation during the 
stopping was less accurate. 
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Figure 7 is a chart, for a fifth exemplary embodiment, based on a 
consist having two locomotives, 50 loaded railcars and 51 empty railcars. Figure 7 
presents the prediction of the same stopping distance using modified initial data. 
Before being used for initialization, the data are adjusted according to the difference 
5 of weight between the two consists. With this new initial data, the prediction is much 
closer to the real stopping distance than in previous determinations. At the end of the 
period, as the model acquires new data, the prediction line approaches the 100% level. 

In Figures 3 through 7, the predictions presented are of the stopping 
distance associated with full brake application that happens at the end of the 2006 
1 0 second period. However, the ATM is capable of providing the stopping distance as 
the brakes are applied. 

Figiire 8 is a chart, for a sixth exemplary embodiment, based on a 
consist having two locomotives, 50 loaded railcars and 51 empty railcars. Figure 8 
presents the predictions of the stopping distance as the brakes are fully applied 
1 5 instantly. In other words, it is the safest distance between two consists following each 
other on the same track. The predictions of the safe distance are made using a model 
that is not initialized. 

Figtire 9 is a chart, for a seventh exemplary embodiment, based on a 
consist having two locomotives, 50 loaded railcars and 51 empty railcars. Figure 9 

20 presents the variation of the predicted stopping distance determined with the ATM 
using initialization information based on data collected after three trips on the same 
track, and using the predictions made on the third trip as references. The difference 
between the safe distance predicted by the ATM that has no initial data and 
predictions made on the third pass of the train on the same track are shown. The 

25 variation of errors for the non-initialized model is between -725 and 1335 feet. 

Figure 10 is a chart, for an eighth exemplary embodiment, based on a 
consist having two locomotives, 50 loaded railcars and 51 empty railcars. Figure 10 
presents the error of the estimated safe distance for the model initialized with 
weighted information. The variation of errors is between -410 and 268 feet. 

30 Model database 

Referring to Figure 1, the data collected by the ATM running on a 
consist is stored in database 16 from where it is retrieved and used as initial data to 
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another consist with the same or close structure. Database 1 6 is accessed on the keys, 
such as, number of locomotives, locomotive type, number of railcars, train length 
(feet), train weight (tons), and load distribution along the consist. 

While the invention has been described in terms of various specific 
embodiments, those skilled in the art will recognize that the invention can be practiced 
with modification within the spirit and scope of the claims. 
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